Velocity correlations and the structure of nonequilibrium hard core fluids 
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A model for the pair distribution function of nonequilibrium hard-core fluids is proposed based 
on a model for the effect of velocity correlations on the structure. Good agreement is found with 
molecular dynamics simulations of granular fluids and of sheared elastic hard spheres. It is argued 
that the incorporation of velocity correlations are crucial to correctly modeling atomic scale structure 
in nonequilibrium fluids. 



The structure of a fluid as characterized by the pair distribution function (pdf) plays a central role in statistical 
mechanics. In equilibrium, it allows one to calculate the equation of state whereas away from equilibrium it is required, 
in the guise of the equal time density-density correlation function, whenever one wishes to project a kinetic equation 
onto the hydrodynamic subspace, see e.g. ref [Q. In both cases, it can be directly measured by light scattering and 
is thus an interesting quantity in order to make connections with experimental studies. The theory of the structure 
of equilibrium fluids is quite advanced and ranges from simple models for spherical hard core systems to perturbative 
models applicable to continuous potentials to integral equations applicable to a wide range of systems ||, [|. Much 
less is known for nonequilibrium fluids. Kinetic theory can yield information in certain regimes such as low density or 
large spatial separations, see. e.g. ||^ but does not provide complete models of nonequilibrium structure applicable 
to dense fluids. Phenomenological models exist based on fluctuating hydrodynamics 0, |^ or Langevin models 
| pj| , Jl2t but these can only be expected to be applicable at large length scales. For moderately dense hard-core fluids, 
the only realistic and tractable kinetic theory is the Enskog equation describing the time-evolution of the one-body 
distribution function (the dense fluid generalization of the Boltzman equation) [nSl . Recently, it has been shown that 
one piece of structural information is directly accessible in hard-core systems within the same set of approximations 
used to derive the Enskog kinetic theory - namely, the pair distribution function at contact . The purpose of this 
Letter is to show that it is possible to use this information to extend the simplest class of equilibrium models, the 
so-called Generalized Mean Spherical Approximation or GMSA Q, to nonequilibrium systems and to validate 
the models constructed against molecular dynamics simulations. In particular, the models will be compared in the 
case of a granular fluid ( a homogeneous system violating energy conservation) as well as strongly sheared elastic 
hard spheres. The latter case is particularly interesting due to the fact that it has motivated previous attempts to 
model nonequilibrium structure, is a realistic model of some sheared colloidal suspensions and has been investigated 
experimentally for such systems [ p^ . 

To construct a model for the structure of a nonequilibrium fluid, we begin, as in nearly all models of equilibrium 
structure, with the Ornstein-Zcrnike (OZ) equation for a uniform fluid 

/i(ri,r2) = c(ri,r2) -f p J drs c{ri,r3)h{r3,r2) (1) 

where p is the density, the structure function /i(ri,r2) is related to the pdf, g(ri,r2), by /i(ri,r2) — 5(ri,r2) — 1 
and eq.(|l]) serves as a definition of the direct correlation function (dcf) c(ri,r2). By definition, hard core atoms 
cannot interpenetrate so the probability of finding two atoms closer together than the hard sphere diameter, cr, is 
zero giving the boundary condition Q (a — ri2) (?(ri, r2) = where ri2 = |ri — r2| and Q{x) is the Heaviside step 
function which is equal to 1 for x > and zero otherwise. Various arguments, based e.g. on the Mayer expansion 
of the various quantities, leads to the conclusion that the direct correlation function is short ranged and the Percus- 
Yevick (PY) theory results from taking it to be zero outside the core: & (ri2 — a) c(ri,r2) — 0. This is sufficient to 
completely specify all quantities and, in particular, the pdf can be analytically determined However, there is 
ambiguity in the equation of state since it may be determined from either the pressure equation, which involves the 
pdf at contact, or from the compressibility equation involving an integral of the structure function over all space: 
these do not yield the same result. A feature of many, if not most, successful theories of equilibrium structure is 
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that they force both routes to the equation of state to yield the same result, a requirement known as thermodynamic 
consistency. The GMS A accomplishes this by giving the direct correlation function a parameterized " tail" by requiring 
O (ri2 — a) c(ri,r2) = u(ri,r2) for some function ■(;(ri,r2). In the original MSA, the tail is taken to be either the 
interatomic potential or the Meyer function constructed from it (thus reducing to the PY approximation in the case 
of hard spheres). In the original GMS A of Waisman jl^, the tail is taken to be a Yukawa w(ri, — Ae~^^^'^'^~^'^ jryi 
and the two parameters are determined by requiring that both the routes to the equation of state yield the Carnahan- 
Starling equation. The pressure equation therefore fixes the pdf at contact whereas the compressibility equation fixes 
its area. This model may be solved analytically and gives a good model of the pdf even for dense fluids and can be 



adapted to give reasonable models other systems such as plasmas |17 . In all cases, the two elements which depend 
on the system being in equilibrium are (1) the arguments for the nature of the direct correlation function outside the 
core and (2) the equations of state and the condition of thermodynamic consistency. 

The first result that suggests the applicability of these models away from equilibrium is the work of Yuste and 



18|. They show that the pdf can be written in terms 
e approximant. This approximant is then taken to be 



Santos on the structure of more complex hard-core systems 
of an auxiliary function which is then parameterized as a Pad( 
the simplest possible that satisfies three physical conditions: the pdf is finite at contact; the pdf goes to one at large 
separations; the Fourier transform of the structure function is finite at all wavevectors. These conditions are sufficient 
to completely specify the auxiliary function and the result is the PY pdf. Extending the approximant by one additional 
term in the numerator and denominator (the relative order of the two is fixed by the constraints) is exactly equivalent 
to the Yukawa closure and the imposition of thermodynamic consistency and the Carnahan-Starling equation of state 
yields the GMSA. By recasting the GMSA solely in terms of general physical requirements and extendable analytic 
approximations, this work accentuates the fact that insofar as the tail of the direct correlation function is simply 
being parameterized, the GMSA is generic. However, before this model can be used in nonequilibrium systems, some 
substitute for the equation of state, used as input, must be found. One role of the equation of state is, as mentioned 
above, to fix the value of the pdf at contact. The second ingredient needed to complete the model is the recent 
result showing that it is possible to calculation any two-body function, evaluated for two atoms in contact, with the 
same degree of approximation as is used in the Enskog equation (the dense fluid generalization of the Boltzmann 
equation for hard core systems) jl^ . A slightly more general form of the result than that given in ref ||l^ is that the 

(nonequilibrium) average value of any two-body function at contact, say J7(r) — {^2ii<_j "^^^3)^ (Qy ~ '^)^ where we 
use the abbreviated notation u(i]) = u(qi, q^, Pi, Pj), is given by 
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where the parameter a characterizes the inelasticity of the hard spheres (a = 1 corresponds to elastic hard spheres), 
f = r/r denotes a unit vector and u'(qi, q^, pi, pj) = ii(qi, q^, p^— ii^q^j (q^ • Py ) , pj-t-ii^qij (q^ • Py )) is the value 
of the two-body function just after the collision has ocured. (This result is easily obtained by applying the method of 
ref. to the pseudo-Liouville equation for inelastic hard spheres ||l^). The notation (...)g indicates an average that it 
is to be performed with the factorized two-body distribution function /2(qj, q^ , Pi, Pj) — Pi)/i(qi, Pj).9o(qj, qj) 

where /i(qi, Pi) is the one-body distribution (normally derived from the Enskog equation) and (7o(qj, q^) is typically 
the local-equilibrium pdf. The first term on the right represents the average over the nonequilibrium state when 
all velocity correlations are neglected whereas the second gives the effect of velocity correlations generated by the 
collision. As discussed previously , the Enskog equation also takes account of these velocity correlations whereas 
both it and eq.(^) neglect correlations that persist for more than one collision. 

In summary, I therefore propose to use a GMSA-like paramterized closure of the OZ equation and to fix the 
parameters using information coming from equation (2) which, in the case u(zj) = 1 gives a generalization of the 
pressure equation applicable to nonequilibrium fluids. The form of the parameterization of the tail of the dcf will be 
guided from results from kinetic theory: the Yukawa form will be used when the pdf is expected to decay exponentially, 
as in the case of a granular system (note that this corrects a previously reported algebraic decay ^ for this system) 
whereas a power- law will be used when the decay is expected to be algebraic (as for shear flow). 

Granular fluids are often modelled as hard spheres which lose energy upon collision in which case the undriven 
system cools at a steady rate. The energy loss per collision, and hence the cooling rate, is controlled by the inelasticity 
parameter, a, introduced in eq.(H) above. In this homogeneous cooling state (HCS) eq.(H) gives g{ri2 — a\a) — ^^^X 
where x is the equilibrium pdf at contact (since the HCS is translationally invariant, all two-body quantities depend 
only on the relative separation). The validity of eq.(||) for this system has been studied in some detail by means of 
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a set of molecular dynamics simulations which will be presented at a later time. The simulations show that it is a 
reasonable approximation in the range 1 > a 0.5 at which point mode-coupling effects become important. Here, we 
only present one comparison, figure (H), between the pdf as determined from simulation and the nonequilibrium GMSA 
with the Yukawa closure for a relatively dense fluid, na^ = 0.5, and a large rate of dissipation, a — 0.5. The model is 
seen to be good not only near the core but also in its representation of the decay of the nonequilibrium contribution 
to the structure. A measure of the significance of the nonequilibrium effects shown is that they are comparable in 
magnitude to the equilibrium pdf which, for this density, varies between a maximum of 2.16 (at the core) and zero, 
with an asymptotic value of 1. In applying the GMSA in this case, we use the known analytic solution of this model 
pot , and must determine the two parameters introduced by the Yukawa tail by the constraints. One constraint comes 
from the value of the pdf at contact but there is no obvious replacement for the compressibility equation used in 
equilibrium so we continue to enforce it, using the expression for the pressure given in ref. 0. Some experimentation 
shows that the model is relatively insensitive to the latter approximation. A more realistic constraint, to be pursued 
in a future publication, would be to replace the compressibility equation by the requirement that the area under the 
pdf agree with that obtained from fluctuating hydrdodynamics calculations [|f^ . 

A second, and more complex test, is the application of the nonequilibrium GMSA to shear flow. Define a cartesian 
coordinate system in which the macroscopic flow field is v(r) = ayx, where a is the shear rate. In this case, an 
expression for the pdf at contact has been studied in some detail [[l4| and, while it is possible to calculate its full 
angular dependence, in the illustrative calculation presented here, we will use the simple approximation g(crr;a) = 
X (1 + ArxTy) with A chosen to give the calculated value at Vx =ry = a/\/2 [Q, which is the dominant contribution. 
To model this, we must use an anisotropic tail for the dcf and therefore take ti(ri, = wo(''i2) + ui(ri2)f 122;? i2j,. The 
solution of the OZ equation with an anisotropic potential is well known , ||^ and will only be sketched here. To begin, 
the structure function and the dcf must be expanded in terms of their angular dependencies: h{ri2) — /iim(?'i2)^im(ri2) 
where repeated indices are summed and Y/„j(r) is a spherical harmonic. We will also need the Fourier transforms 

of these expansions which takes the form /i(k) = himik)Yim(k) where him{k) may be calculated from /i;m(?'i2) via 
a Hankel transformation (in the case I ~ this reduces to a Fourier transform). Fourier transforming the OZ 
equation and multiplying by Y*^ (k) and integrating then gives an infinite set of coupled equations. The minimal set 
of components that we can work with are the spherically symmetric component, /ioo('"i2)j and the components that 
combine to produce the angular dependence of the tail, VyixXviyi namely hi^iTvi) and h^-iirvi)- Furthermore, we 
know that at contact, these give the dominant contributions to hivyi). For illustrative purposes, we therefore keep 
only these three components and the corresponding equations in the hierarchy with the result that we must solve the 
system 

/loo = Coo + '^\f^ [coo/Joo + C22/12-2 + C2-2/122] 

/i22 = C22 + nJ — (cooft.22 + C22/100) 
V 47r ' 

/12-2 = C2-2 + "\/ T- (coo/i2-2 + C2-2/100) (3) 

together with the boundary conditions described above. Furthermore, a consideration of the symmetry of the problem 
and the form of the boundary conditions shows that ft-22 must be pure imaginary and that ^,22 — ^2-2 that the 
independent equations can be written as 

(/loo ± »/l22) = (coO ± jC22) + ^^^^ i^OO ± i/122) (coO ± iC22) (4) 

which has the form of two decoupled OZ equations (at this point, the solution maps to that of the dipolar hard-sphere 

model [|). The fuU pdf now has the form ^(r) = 1 -f- --^/loo(r) + i^J^h22{'r)YxYy. The only element that remains 

is to specify the tails of the dcf. At contact, we know that /ioo('') should be unchanged from its equilibrium value 
whereas h22{r) is substantial. This suggests that we take vo{ri2) = although in a more sophisticated treatment 
we might want to use a Yukawa so as to fix the value of /loo(^) at contact (from the model for the pdf at conact, 
this should be unchanged from equilibrium but eq.(^ will not necessarily respect this - however, upon solving the 
model, the change of this component at contact is found to be minimal). For the anisotropic part, the simplest 
choice we can make is to take fi(ri2) ~ B/rf2 (since then these equations may be solved using the FY results [^). 
Although kinetic theory calculations H] indicate the presence of a long-ranged tail for h22{r) that decays as l/r, we 
nevertheless proceed using vi (ri2 ) = B /r J2 because (a) even using a longer-ranged I /r tail for the dcf will not yield 
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1/r behavior in the pdf |17| and (b) no matter what tail we choose, there will be a decay The constant B 
is determined so as to give agreement with the predicted value of the pdf at contact and since there are no other 
parameters, the ad hoc use of the compressibility equation is not needed in this case. The solution then follows closely 
that for dipolar hard spheres, see e.g. To test this, I have performed molecular dynamics computer simulations, 
following the procedures described in ref . jl^ , to determine the functions /iqq (r) and /122 (f) and a comparison to the 
nonequilibrium GMSA model is shown in Figs. (2) and (3) for a moderately dense fluid, na'^ = 0.5, in a strongly 
driven state, aioo — 0.5 where too — cr/ {An* \/TrkBT) is the Boltzmann collision time. For hoo{r), there is little 
difference between the theoretical result, which is virtually identical to the FY pdf, and the simulation indicating 
little or no change to this projection of the pdf due to the shear flow. For the more interesting h22{r) component, 
the agreement is not exact but this crude calculation nevertheless gives a reasonable approximation to the amplitude 
and peak of the oscillations of this quantity which arises solely from the nonequilibrium state and which are seen by 
comparison to fig. (2) to be comparable in magnitude to the equilibrium pdf. A more complete solution to this model 
is in preparation. 

These results on two very different strongly nonequilibrium systems show that the knowledge of the velocity correla- 
tions produced when two atoms collide, and the consequent contribution to the pdf, is sufficient to generalize a simple 
class of structural models to nonequilibrium fluids. While these models may be useful as input to kinetic theory or 
for the interpretation of light-scattering experiments, the most important point to be made here is the critical role 
played by the velocity correlations in determining the atomic-scale structure in nonequilibrium systems. Theories 
of nonequilibrium structure which do not explicitly take into account atomic-scale velocity correlations, such as the 
phenomenological theories of Hess and Ronis and the semi-phenomenological theory of Eu |Q, are, at 

least when applied to hard-core systems, a priori unlikely to be able to model the small-length scale distortions of 
the fluid structure produced by nonequilibrium processes while, as shown above, even the simplest theories which do 
take them into account compare reasonably well to simulation. 
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Figl 




FIG. 1. The difference between the pdf for a model granular fluid with na^ = 0.5 and a = 0.5 and the PY pdf for an 
equilibrium fluid of the same density. The circles are from simulation and the line is the GMSA result. 
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Fig. 2 
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FIG. 2. The quantity feoo('") = ':^^oo{r) for a fluid of sheared hard spheres with density na^ = 0.5 and reduced shear rate 
atoO = 0.5 where toO is the Boltzmann colhsion time. The circles are from simulation and the line is the GMSA result. 
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Fig. 3 




FIG. 3. Same as Fig. 2 except showing the quantity /122 = i\/^h22{r). 
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